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Abstract 



^D \ A punctuated equilibrium model of biological evolution with relative fitness between 

^^ ' different species being the fundamental driving force of evolution is introduced. Mu- 

tation is modeled as a fitness updating cellular automaton process where the change 

^" ■ in fitness after mutation follows a Gaussian distribution with mean re > and 

standard deviation a. Scaling behaviors are observed in our numerical simulation, 
J^ \ indicating that the model is self-organized critical. Besides, the numerical experi- 

ment suggests that models with different x and a belong to the same universality 

'^ ' class. 

a" 

o 

o 

> 

X 



FAGS numbers: 87.10. -Fe, 05.40.-hj 



Typeset using REVI^ 



*Fresent Address: School of Natural Sciences, Institute for Advanced Study, Olden Lane, Frinceton, NJ 
08540. 



Biological evolution takes place in bursts separated by relatively long periods of quiescence ^ 
In fact, extinction may be episodic at all scales 0]. To capture the essence of this scaling behavior, 
a punctuated equilibrium model of biological evolution is introduced recently by Bak and Sneppen 
(BS Model) 0. An ecosystem in the model is made up of N species, each associated with a scalar 
between zero and one, called fitness of the species. The higher the fitness, the more adaptive to the 
ecosystem and hence the less likely to mutate the species is. Thus, the species with the minimum 
fitness in the system is most likely to mutate or extinct next. Fitness is therefore also a measure of 
the barrier against mutation. 

In the one-dimensional BS model, the A^ species are arranged on a line with periodic boundary 
condition. At each timestep, the species with minimum fitness and its two nearest neighbors 
mutate; their fitnesses are replaced by uniform and uncorrelated random numbers between zero and 
one 0. The updating process is repeated forever. Scaling behavior is observed in the distribution 
of distance between successive mutations once the stationary state is reached, which is a signature 
of self-organized criticality. Generalization to the A^- dimensional and the mean field models can be 
made, and the results are essentially the same as their one-dimensional counterpart [|,0|. 

However, there are a number of conceptual difficulties in the BS model. First, it is unclear why 
mutation changes the fitness of the neighboring species and why such changes take place all at the 
same time when the least-fit species mutates. This kind of "action at a distance" couplings between 
the least-fit species and its neighbors are not physically justified. Second, there is no correlation 
between the fitnesses of a species before and after mutation (or extinction). This is not consistent 
with the fundamental assumption of the Darwinian school that all existing species share a common 
ancestor p . Finally, once the stationary state is reached, the average fitness of the species becomes 
a constant. Thus species in the present time and those in the past have the same average adaptive 
power to the environment, which is not true in nature. 

Now we are going to introduce an alternative model of biological evolution free from all the 
above conceptual difficulties. Our ecosystem is made up of a single food chain consisting of A^ 
species labeled from 1 to A^, with 1 being the producer, and A^ being the terminal consumer. The 
fitness of the i-th species is denoted by /«. Similar to the case of energy in physics, the value of 
fitness is meaningful only when it is measured with respect to a reference fitness level. The actual 
fitness value of a species is not a physical observable because the value of the reference fitness can be 
arbitrarily assigned. The relative fitnesses between different species (together with the environment) 
are the only physical measurable and meaningful quantities when describing how much adaptive 
power the species have. 

In general, the "relative fitness" of a species can be written as 

V^ = m^-fE)+T.Mf^-fk) (1) 

for some at-, (3k > where /g is the "fitness" of the environment. The first and second terms in 
Eq. (|l|) represent the competition of the species with the environment and with the other species 
respectively. In our present model, only the special case where Vi = 2/j — /j+i — /j_i is considered. 
Our "nearest neighbor interaction assumption" is justified provided that competition is only im- 
portant between the predator and the prey. Further analysis in the case where the environmental 
factors cannot be neglected is underway 0. 

Relative fitness can be also regarded as a measure of the barrier against mutation [Q. Provided 
that the probability of mutation of a species in unit time is history independent, the characteristic 



mutation time for this species is given by: 

g{v) = Ae''' (2) 



for some A, \ > 0. The above interpretation is consistent with the NK Model and its variations |]T0 
in which mutation is regarded as hopping between nearby locally optimal "fitness landscape". The 
characteristic time required for mutation therefore depends exponentially on the activation barrier. 
Given the fact that some species mutate in a much faster time scale than others, we expect A ^ 
PJTT[]. With a sufficiently large A, the species with the smallest relative fitness is almost surely to 
be the next one to mutate or extinct. A lower bound for A will be given later in the text. 

With this assumption in mind, the species i with minimum relative fitness is always chosen to 
mutate by updating its fitness fi as: 

where B^^j is a random number chosen from a Gaussian distribution with mean x and standard 
deviation a. We choose x > to refiect the trend that newly mutated species is likely to be 
more adaptable to the ecosystem. Moreover, we require a ^ a; so as to allow exceptions to the 
general trend. The relative fitnesses of all species in the ecosystem are then recalculated. Clearly, 
only the relative fitnesses of species i — 1, i, and i + 1 will change after the mutation of i. A 
single cellular automaton timestep is now completed, and the process is repeated indefinitely. The 
above construction is possible if we assume that the mutation (or replacement) of a species in our 
ecosystem does not alter the structure of the food chain. Unlike the BS model, correlation between 
new and old fitnesses of the mutating species is properly accounted for. In our model, the change in 
evolutionary landscape of a neighboring species j is a result of the coupling between the mutating 
species i and its neighbor j via the change in relative fitness Vj. The absolute fitness /,■ of the 
species remains unchanged right at the time of the mutation oi i. In this way, relative fitnesses may 
be regarded as "fields" which couple between neighboring species, telling them the status of their 
neighbors. Our "particle-field-particle" interaction model looks more natural than the "action at a 
distance" BS model. 

We complete the model using the closed boundary conditions, namely, vi = /i— /2 and vn = /at — 
/n-i- Closed boundary conditions are consistent with the fact that /, are not physical observables: 
once we know the values of fj, fi can only be determined up to a common additive constant. Finally, 
we would like to emphasize that the real physical time elapsed in each cellular automaton timestep 
is not a constant and depends on the minimum relative fitness of the system P,|ll|] . 



By working with the relative fitness instead of the absolute one, Eq. (0) can be rewritten as 



'new 



i^i~l)old - Bx,a 



{'"i)new = (^Oo;d + 25.,. (4) 



■'X,(7 



except possibly at the system boundaries. Eq. (^) is similar to the updating rules used by Bak and 
Sneppen in [^, with the important exception that the random variables we used for all the three 
species are correlated. We therefore expect our model to exhibit self-organized critical phenomena 



because it is observed in the "less correlated" BS model. This is indeed the case as we proceed to 
show it numerically. 

Numerical simulation is carried out with A^ = 1024. We fix x to be 1.0, and systems with other 
values of x can be rescaled to ours by a simple scale transformation. We have considered cases 
with a between 0.1 and 3.0. In all cases, after a long transient period, the system settles into a 
stationary state independent of the initial fitnesses of the species. Then a further 2 x 10® mutations 
are recorded for statistical purpose. The distribution of the relative fitness v and the minimum 
relative fitness Vmin are shown in Fig. |l] and Fig. ^ respectively. Similar to the BS model, there 
is an upper critical value Vc above which there is zero probability of the minimum relative fitness 
Vmin- Although the precise values of Vc (see Table D is a dependent, the shapes of the distribution 
curves do not change. Unlike the BS model, the distribution of v (and Vmin) does not fiatten when 
its value is greater than (or less than) Vc- In fact. Fig. 0b tells us that (except for the small hump 
near Vc) the distribution of Vmin {D{vmin)) can be well approximated by 

in) l-t '^min ^ '-'c /|-\ 

otherwise 

for some B, u > 0. From Eqs. (^ and (^), the distribution of time between successive mutations 
T>{T) is then given by [O]: 




X,(r)^^Z,(i,4)«f(i)Vi- (6) 

whenever T < Ae^^". 

The average mutation rate of a species is J D{v)/g{v)dv, and the minimum possible mutation 
rate of the least fit species is l/g{vc)- The assumption of always choosing the least adaptable species 
to mutate is justified provided that the sum of mutation rates of all the other species is much less 
than the mutation rate of the least adaptable one |TT[; thus (A^ — 1) / D{v)/g{v)dv <^ l/g{yc). This 



is true if A ^ A^. So for the case where A^ is sufficiently large, z//A is negligible and hence 1/T 
scaling is expected for V{T). 

Suppose i and j are the mutating species at the present and the previous timesteps respectively. 
We define the distance from the last event d as |^— j|. The distribution of d {F{d)) is shown in Fig. ^ 
where l/rf" scaling is observed for about 2.5 decades with a = 2.87 ± 0.02, which is independent 
of the values of x and a used in all of the simulations. So we strongly suspect that the critical 
exponent a is independent of x and a, and models with different x and a in fact belong to the 
same universality class. Furthermore, scaling behavior in both V{T) and F{d) are clear signs of 
self-organized criticality, that is, under its own dynamics, the ecosystem is driven into a state where 
there is no characteristic correlation length and time scales. 

For a given value of relative fitness Vref, an avalanche of size s is defined to be the process 
of s consecutive mutations with minimum relative fitnesses less than or equal to Vref, which are 
immediately preceded and followed by mutations with minimum relative fitnesses greater than Vref 
^. For sufficiently large A, the distribution of avalanche size coincides with that of the total number 
of mutations occurred in a time interval of g{vref) Jill- Fig- § shows that 1/s^ scaling holds over 



three decades; the value of the critical exponent /? is found to be 0.99±0.02 and independent of the 
value of a. 

In summary, we have introduced a self-organized critical model of biological evolution in which 
competition between between species is the driving force using the idea of relative fitness. Relative 



fitness acts like a background field and tells the species to response to its own situation relative 
to the other species and to the environment. Although the average value of relative fitness is a 
constant in the system, the average absolute fitness increases with time, reflecting the improvement 
of adaptive ability. The time-average relative fitness of a species is constant which is consistent 



with the Red-Queen hypothesis |[T2|. Our numerical experiments suggest that models with different 
values of x and a belong to the same universality class. Our model may also belongs to the same 
universality class as the Robin Hood model whose critical exponent a = 2.86 ±0.01 |T^. Numerical 
experiments using periodic boundary conditions (fi = 2/i — /2 — /at, vjy = 2/jv — /i — Jn-i), and 
reflective boundary conditions {vi = 2/i — 2/2, vjy = 2/jv — 2fiy_i) are also performed. We found 
that the shape of all the distribution curves, together with the critical exponents, are insensitive to 
the boundary conditions used. 
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TABLES 

TABLE L Upper critical relative fitness Vc as a function of a. In all the cases listed below, N = 1024 
and X = 1.0. 

o^ Vc 

0.5 -0.44±0.01 

1.5 -0.72±0.01 

3.0 -1.23±0.01 



FIGURES 

FIG. 1. Distribution of relative fitness v. a = 0.5, 1.5, and 3.0 are used in the dotted, dash, and solid 
curves respectively. This convention is used in all subsequent figures. 

FIG. 2. Distribution of minimum relative fitness Vmin- This is shown in linear plot in (a) and in semilog 
plot in (b). 

FIG. 3. Distribution of the distance from last event d. Power law behavior is observed over 2.5 decades 
in all the three cases, and their critical exponents agree with each other and are equal to 2.87±0.02. 

FIG. 4. Distribution of avalanche size below a reference fitness level Vref, which is chosen in the way 
that about 98% of the minimum fitnesses will fall below it. These values correspond to —1.4, —0.8 and 
—0.47 for the solid, dash, and dotted curves respectively. Scaling behavior is observed over 3 decades in 
all the three curves. The scaling region will lengthen as v^ef approaches Vc- The critical exponent /? is 
measured to be 0.99±0.02 in the regions where the power law holds. 



